# Get Table 6
# =========================================================================================================
# remove everything
rm(list=ls())

# =================================================================================

# Load packages
library(foreign)

# Load the bootstrapped data
# To get the full tables 6 and I5, names of the following two files need to be changed for each model: F1, F2, F3, F4, W1, P1, P2, P3, IPW, IV, Lewbel
load("Boot_s_Dt_F1.RData")
load("Jack_s_Dt_F1.RData")

# ---------------------------------------------------------------------
# Get the bootstrapped CI.
B <- 400

# ================================================================================================
# ================================================================================================

# summary of the estimates
SPE <- summary(Bpsw[,1])
# summary of the Bootstrapped estimates
SPB <- apply(Bpsw, MARGIN = 2, summary)

# acceleration factor
SPE.jk <- apply(JKpsw, MARGIN = 2, summary)
acc.top <- apply((apply(SPE.jk, MARGIN = 1, mean) - SPE.jk)^3, MARGIN = 1, sum)
acc.bot <- 6*(apply((apply(SPE.jk, MARGIN = 1, mean) - SPE.jk)^2, MARGIN = 1, sum))^(3/2)
acc.factor <- acc.top/acc.bot

# ---------------------------------------------------------------------
# Median: the 3rd position of SPE is the mean estimate
pos <- 3
dec <- 4
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.md <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# Percentage of positive SP obs.
dec <- 4
phio <- qnorm(apply((Bpsw[,1:400] < Bpsw[,1]), MARGIN = 1, sum)/B)

phi5 <- qnorm(0.05)

acc.top2 <- apply((apply(JKpsw, MARGIN = 1, mean) - JKpsw)^3, MARGIN = 1, sum)
acc.bot2 <- 6*(apply((apply(JKpsw, MARGIN = 1, mean) - JKpsw)^2, MARGIN = 1, sum))^(3/2)
acc.factor2 <- acc.top2/acc.bot2
c <- acc.factor2
a12 <- pnorm(phio + (phio+phi5)/(1-c*(phio+phi5)))
up <- c()
for(jj in 1:nrow(Bpsw)){
  temp <- as.numeric(quantile(Bpsw[jj,1:400], probs = a12[jj]))
  up <- rbind(up, temp)
  if(jj==round(jj,-3)) print(jj)
}
per <- format(round(sum(up>0)/nrow(Bpsw), dec), nsmall = dec)

# ---------------------------------------------------------------------
# Tabulate
print(paste(SPE[3],per, sep = ";"), quote = F)
print(paste(ci.md, sep = ";"),quote = F)
sp <- median(Bpsw[,1]); per1 <- per; ci.md1 <- ci.md


# ================================================================================================
# ================================================================================================

# summary of the estimates
DLE <- summary(BpF[,1])
# summary of the Bootstrapped estimates
DLB <- apply(BpF, MARGIN = 2, summary)

# acceleration factor
DLE.jk <- apply(JKpF, MARGIN = 2, summary)
acc.top <- apply((apply(DLE.jk, MARGIN = 1, mean) - DLE.jk)^3, MARGIN = 1, sum)
acc.bot <- 6*(apply((apply(DLE.jk, MARGIN = 1, mean) - DLE.jk)^2, MARGIN = 1, sum))^(3/2)
acc.factor <- acc.top/acc.bot

# ---------------------------------------------------------------------
# Median: the 3rd position of DLE is the mean estimate
pos <- 3
dec <- 4
phio <- qnorm(sum(DLB[pos,] < DLE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(DLB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(DLB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.md <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# Percentage of positive DL obs.
dec <- 4
phio <- qnorm(apply((BpF[,1:400] < BpF[,1]), MARGIN = 1, sum)/B)

phi5 <- qnorm(0.05)

acc.top2 <- apply((apply(JKpF, MARGIN = 1, mean) - JKpF)^3, MARGIN = 1, sum)
acc.bot2 <- 6*(apply((apply(JKpF, MARGIN = 1, mean) - JKpF)^2, MARGIN = 1, sum))^(3/2)
acc.factor2 <- acc.top2/acc.bot2
c <- acc.factor2
a12 <- pnorm(phio + (phio+phi5)/(1-c*(phio+phi5)))
up <- c()
for(jj in 1:nrow(BpF)){
  temp <- as.numeric(quantile(BpF[jj,1:400], probs = a12[jj]))
  up <- rbind(up, temp)
  if(jj==round(jj,-3)) print(jj)
}
per <- format(round(sum(up>0)/nrow(BpF), dec), nsmall = dec)
print(per, quote = F)

# ---------------------------------------------------------------------
# Tabulate
print(paste(ci.md, sep = ";"),quote = F)
dl <- median(BpF[,1]); per2 <- per; ci.md2 <- ci.md

# ================================================================================================
# ================================================================================================

Bptil <- Bptil[complete.cases(Bptil),]; JKptil <- JKptil[complete.cases(JKptil),]
# summary of the estimates
TILE <- summary(Bptil[,1])
# summary of the Bootstrapped estimates
TILB <- apply(Bptil, MARGIN = 2, summary)

# acceleration factor
TILE.jk <- apply(JKptil, MARGIN = 2, summary)
acc.top <- apply((apply(TILE.jk, MARGIN = 1, mean) - TILE.jk)^3, MARGIN = 1, sum)
acc.bot <- 6*(apply((apply(TILE.jk, MARGIN = 1, mean) - TILE.jk)^2, MARGIN = 1, sum))^(3/2)
acc.factor <- acc.top/acc.bot

# ---------------------------------------------------------------------
# Median: the 3rd position of TILE is the mean estimate
pos <- 3
dec <- 4
phio <- qnorm(sum(TILB[pos,] < TILE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(TILB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(TILB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.md <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# Percentage of positive DL obs.
dec <- 4
phio <- qnorm(apply((Bptil[,1:400] < Bptil[,1]), MARGIN = 1, sum)/B)

phi5 <- qnorm(0.05)

acc.top2 <- apply((apply(JKptil, MARGIN = 1, mean) - JKptil)^3, MARGIN = 1, sum)
acc.bot2 <- 6*(apply((apply(JKptil, MARGIN = 1, mean) - JKptil)^2, MARGIN = 1, sum))^(3/2)
acc.factor2 <- acc.top2/acc.bot2
c <- acc.factor2
a12 <- pnorm(phio + (phio+phi5)/(1-c*(phio+phi5)))
up <- c()
for(jj in 1:nrow(Bptil)){
  temp <- as.numeric(quantile(Bptil[jj,1:400], probs = a12[jj]))
  up <- rbind(up, temp)
  if(jj==round(jj,-3)) print(jj)
}
per <- format(round(sum(up>0, na.rm = TRUE)/nrow(Bptil), dec), nsmall = dec)
print(per, quote = F)

# ---------------------------------------------------------------------
# Tabulate
print(paste(ci.md, sep = ";"),quote = F)

til <- median(Bptil[,1]); per3 <- per; ci.md3 <- ci.md

format(round(sp, dec), nsmall = dec); print(ci.md1, quote = FALSE); print(per1, quote = FALSE)
format(round(dl, dec), nsmall = dec); print(ci.md2, quote = FALSE); print(per2, quote = FALSE)
format(round(til, dec), nsmall = dec); print(ci.md3, quote = FALSE); print(per3, quote = FALSE)


